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We present a methodology for modelling population dynamics with formal means of computer sci- 
ence. This allows unambiguous description of systems and application of analysis tools such as 
simulators and model checkers. In particular, the dynamics of a population of Aedes albopictus (a 
species of mosquito) and its modelling with the Stochastic Calculus of Looping Sequences (Stochas- 
tic CLS) are considered. The use of Stochastic CLS to model population dynamics requires an 
extension which allows environmental events (such as changes in the temperature and rainfalls) to 
be taken into account. A simulator for the constructed model is developed via translation into the 
specification language Maude, and used to compare the dynamics obtained from the model with real 
data. 

1 Introduction 

In the last few years many formalisms have been defined to model biological systems at molecular 
and cellular levels 0|9l[T3]|22l|23]]. These formalisms allow unambiguous description of systems and 
application of analysis tools, such as simulators and model checkers. 

Among these formalisms the Calculus of Looping Sequences (CLS) Q seems to be applicable to 
other classes of biological systems. CLS is based on term rewriting, in which terms may represent simple 
biological structures and compartments, and rewrite rules may represent very general events. Moreover, 
a stochastic extension of CLS has been defined, called Stochastic CLS [2], which allows the dynamics 
over time of the described system to be studied (6j|7l. 

In this paper we deal with the problem of modelling population dynamics with formal means of 
computer science. Many aspects of population dynamics, such as births, deaths and interaction of indi- 
viduals, can be modelled by using Stochastic CLS. Other aspects related to environmental events, such 
as changes in climatic conditions, require an extension of the formalism. In this paper we define such an 
extension and use it to model the dynamics of a population of Aedes albopictus. 

Aedes albopictus (Skuse), or Asian tiger mosquito, is a species indigenous to the oriental region, but 
it is now widespread in many countries throughout the world. It is an aggressive mosquito, which causes 
nuisance and it is well known as an important disease vector fl9l . 

P. Milazzo and M.J. Perez Jimenez (Eds.): Applications of Membrane Computing, © T.A. Basuki et al. 
Concurrency and Agent-based Modelling in Population Biology (AMCA-POP 2010) This work is licensed under the 
EPTCS 33, 2010, pp. 184361 doi ll0.4204/EPTCS.33.2l ICreative CommonslAttributionl License. 



T.A. Basuki et al. 



19 



A simulator for the constructed model is developed via translation into the specification language 
Maude [12], and used to compare the dynamics obtained from the model with real data. 

There are a number of other approaches to the modelling of population dynamics with formal means 
of computer science lISTl 1011201 . Barbuti et al. Q extend P systems with features typical of timed automata 
with the aim of describing periodic environmental events such as changes of seasons. Cardona et al. iflOl 
propose a modelling framework based on P systems and apply it to the modelling of the dynamics of 
some scavenger birds in the Pyrenees. McCaig, Norman and Shankland [20] present a process algebraic 
approach to the modelling of population dynamics. With respect to these proposals we believe that our 
approach allows a finer modelling of environmental events. Moreover, thanks to the extensions of the 
tools already developed for Stochastic CLS, it offers means for accurate analysis of phenomena. 

2 Stochastic CLS 

Calculi of Looping Sequences (CLS class) is a class of formalisms introduced in Milazzo's PhD the- 
sis 11211 for modelling biological systems. The first formalism of the class to be defined, the Full Calcu- 
lus of Looping Sequences (Full CLS) uses 4 operators: sequencing, parallel composition, looping and 
containment. The parallel composition operator has the typical semantics as in other formalisms such as 
the 71-Calculus and Brane Calculi. Sequencing is inspired by the sequential structure of several macro- 
molecules such as DNA. The looping operator is always applied together with the containment operator 
and supports the modelling of membrane-like structures. An important language of the CLS class is 
Stochastic CLS, which supports the modelling of quantitative aspects of biological systems such as time 
and probabilities. 

We start by introducing the syntax of sequences and terms, the basic building blocks of Stochastic 
CLS. 

Definition 1. Sequences S and Terms T are defined as follows: 

S ::= £ | S-S | a T ::= S | (T) L J T \ T \ T 

where £ represents the empty sequence and a^<§. We denote the set of all terms with , and the set of 
all sequences with .9 1 . 

We assume the existence of a possibly infinite set of symbols <§ . The parallel composition operator 
| is used to model a mixture of elements. The application of the looping and containment operator to 
two terms T\ and T2, denoted by (Ti) L j T2, models structure T2 within a compartment surrounded by 
structure 7\ . Structure T\ is called the loop part and T2 is called its content part. 

The behaviour of a biological system is modeled by means of transitions between terms. This is done 
by applying rewrite rules, that are described by two patterns, to be instantiated by terms, and a rate that 
defines the frequency with which the rule is applied. 

Definition 2. Let TV be an infinite set of term variables ranged over by X ,Y,Z, . . . Term Patterns TP 
and Patterns P are defined as follows: 

TP ::= 5 | (P) L J P \ TP\TP P ::= TP \ TP\X 

where X € TV. We denote with & the set of all patterns. We denote with Var(P) the set of variables in 
P. 
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Definition 3. An instantiation is a partial function o : TV — > 2f. We denote with £ the set of all possible 
instantiations. Given P G we denote with Po~ the term obtained by replacing all variables X € Var(P) 
with o(X). 

Definition 4. A rewrite rule is a triple (Pl,Ic,Pr), denoted with Pl^Pr, where Pl,Pr € £P, k € R and 
such that Var(P R ) C Var(P L ). 

Definition 5. A biological system is a pair {T,M), where T is a term representing the initial state of the 
system and 2% is a set of rewrite rules that represent the potential events which may occur in the system. 

Interactions between populations in a biological or ecological system occur through some kind of 
reactions, which may be biochemical reactions at molecular level or changes in organisms' development 
in ecosystems. To perform in silico analysis of a biological or ecological system, the behaviour of the 
system must be simulated (in silico experiment). The problem of simulating chemically reacting system 
was stated by Gillespie [17]. We generalise Gillespie's formulation of the problem to any biological or 
ecological system as follows. 

A volume or environment V contains a mixture of N species S\,... ,Sn which can inter- 
react through M reaction channels (R\,. . . ,Rm)- Given the initial numbers of individuals 
(molecules or organisms) of each species, what will these population levels be at any later 
time? 

Gillespie consider time evolution of a reacting system as discrete and stochastic. In Gillespie's Stochastic 
Simulation Algorithm ifTol . the state of the system is represented by a vector x = X(?) = (X\ (t),--- ,XN(t)), 
where Xj(t) represents the number of Sj individuals in V at time t. Gillespie assumed that for every 
reaction channel Rj, there is a constant cj such that cjdt is the average probability that a particular com- 
bination of reactant individuals in Rj will react accordingly in the next infinitesimal time interval dt. To 
calculate the probability that a reaction Rj will occur in V in the next infinitesimal time interval (t , t +dt), 
we must multiply Cjdt by the total number of distinct combinations of individuals in V at time t that are 
reactants of Rj. Let us denote such number by hj(x). Gillespie defines the propensity function a,(x) for 
reaction Rj as the product of hj(x) and cj, such that aj(x) dt is the probability that one Rj reaction will 
occur in the next infinitesimal time interval [t,t + dt). 

Gillespie defines a Direct Method to implement his Stochastic Simulation Algorithm [ 171 . This 
version of Gillespie's SSA is defined as follows. 

Algorithm 1. Let {R\ , . . . ,Rm} be a set of rewrite rules, X\,...,Xn be numbers of ,/V categories of 
individuals, maxtime be the time limit for the duration of the simulation. 

Step Initialise simulation time t to 0. Compute propensity a, for every rewrite rule Rj. 

Step 1 Compute the time increment T. 

Step 2 Increase simulation time t by time increment T. 

Step 3 If t > maxtime then stop. Otherwise select the next rule index /J.. 

Step 4 Execute rule R^ and update numbers X\,. . . ,X^ of N categories of individuals and propensities 
ai for all rewrite rules Rj affected by the application of R^ accordingly. Return to Step 1. 

Gillespie showed in his paper |fT7l that the time when next reaction occurs (time t + t calculated at 
Step 2 when reaction R^ selected at Step 3 occurs) is exponentially distributed with parameter «o( x ) = 
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Ljti o-i ( x )- Gillespie used a general Monte Carlo method called inversion method to compute the expo- 
nentially distributed X and jj. from two uniformly distributed random numbers as follows: 



where r\^r% are two real values uniformly distributed over interval [0,1] generated by a random number 
generator. 

In previous work Basuki, Cerone and Carvalho [6] extended Algorithm Q] to handle compartment 
selection. This is useful when we have to simulate a biological system with multi-compartments as is 
the case for molecular reactions occurring within cells. Such a modified version of the Direct Method is 
described as follows. 

Algorithm 2. Let {R\,. . . ,Rm} be a set of M rewrite rules, X\,...,Xn be numbers of N categories of 
individuals, maxtime be the time limit for the duration of the simulation. 

Step Initialise simulation time t to 0. Compute propensity a, for every rewrite rule Rj. 

Step 1 Compute the time increment T. 

Step 2 Increase simulation time t by time increment T. 

Step 3 If t > maxtime then stop. Otherwise select the next rule index pt and the index 6 of the compart- 
ment in which rule R^ will occur. 

Step 4 Execute rule R^ in the compartment with index 6 and update numbers X\,... ,X^ of N cate- 
gories of individuals and propensities a, for all rewrite rules /?, affected by the application of R^ 
accordingly. Return to Step 1. 

Since reactions are confined within compartments, we need to extend Gillespie's algorithm to choose 
in which compartment reaction R^ should occur. Let C be the number of compartments and X' k the 
number of individuals of kind St in the z-th compartment. We define = J^LjXL 

Let a'j be the propensity of reaction Rj occurring inside the z'-th compartment. Then a'- is defined as 
the product of cj by the number h 1 - of distinct combinations of reacting individuals of reaction Rj within 
the z'-th compartment. We define aj = Y%=\ a )- If t i s tne current simulaton time, then t + r represents the 
time at which next reaction occurs, with % exponentially distributed with parameter ao = Y?f=\ a j- Time 
increment % is calculated as in Alghorithm Q] The index \i of the reaction that occurs at time t + % and 
the index 6 of the compartment in which such reaction occurs are calculated as follows: 




(1) 



H-l n 
}JL = the integer for which V a v (x) < r2ao(x) < V a v (x) 



(2) 




m o-i /j e 

(lJ.,d) = the integers for which V a) < r2ao < }, / , a) 

;=li=l j=li=l 



(3) 



where ri is a random real number which is uniformly distributed over interval [0,1]. 



3 Extending Stochastic CLS 



The evolution of a system modelled by using Stochastic CLS is entirely characterised by the rewrite 
rules, which determine the occurrence of events in the system. In this way the set of rewrite rules 
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predicts all events that may occur. This works well for biological systems, where all events are caused 
by biochemical reactions which are governed by precise laws. 

In ecological systems, instead, we need to deal with environmental events, whose cause is often un- 
known or depends on a very complex combination of factors, which are external to the system itself. For 
example the dynamics of a population of a given species depends not only on the interaction with other 
species within the same ecosystem, such as predators, preys and competitors, but also on the occurrence 
of environmental events such as climatic events (i.e. variation of temperature and rainfalls) and events 
related to habitats (i.e. tree clearing, desiccation of a water container, pollution, hunting and human set- 
tlement). Therefore, we assume the existence of a list of external events, with information about the 
time when these events occur. The occurrence of an external event may modify some environmental 
information which affects the ecosystem evolution, such as temperature, volume of water, desiccation, 
level of pollution. Moreover, the list of external events may change dynamically. For instance, an initial 
desiccation event for a water container will be removed from the list after the occurrence of a rainfall 
event, and will be replaced with a new desiccation event with a later desiccation time. 

We extend Stochastic CLS by introducing a list E of external events. The events in list E are sorted in 
increasing order based on the time they are scheduled to occur. When an external event occurs it updates 
information in the system state. The updated information may be then used by rewrite rules. 

In general, the environment is organised as several nested compartments, each associated with spe- 
cific environmental information, which is relevant to the specific ecosystem we are modelling and may 
be modified by the occurrence of external events. We further extend Stochastic CLS by attaching en- 
vironmental information to the looping operator. This is similar to the extension of Stochastic CLS to 
Spatial CLS |4), in which spatial information is added to the looping operator and sequence. 

Definition 6. Terms T, Nonparallel Terms C, Sequences S and Environmental Information / are given 
by the following grammar: 

T ::= C" | T\T C ::= S \ (r)J'J T 

S ::= £ | a | S-S I ::= X \ a:V | // 

where a is a generic element of <§, £ represents the empty sequence, A represents the empty environmental 
information, V represents the information value and ngR We denote with , c €, 5? and the infinite 
set of terms, nonparallel terms, sequences and environmental information, respectively. 

Note that in Definition [6] we have introduced a notation to group identical nonparallel terms together. 
For instance, C 5 is equivalent to C \ C \ C \ C \ C. 

Events in event list E update environmental information in the system state. Every element of E is 
a triple (Ne^Ve^e), where Ne is the name of the event, Vg is a value that will be used to update the 
information field related to this event and ?£ is the time at which this event is scheduled to occur. We 
assume the existence of an event handler algorithm which will handle the update of the term representing 
the system state due to the occurrence of an event (Ne, Ve^e)- 

To run the simulation using the extended version of Stochastic CLS, we need to modify the version 
of Gillespie's Direct Method [ 17] defined in Algorithm |2l In modelling population dynamics we have to 
deal with the same problem we encounter at cellular level: reactions occur in compartments. Therefore, 
we extend the Direct Method for multi-compartments described in Algorithm [2] with additional steps to 
handle the execution of external events from the event list. After computing the time of the next rewrite 
rule, we need to compare this time with the time of the first event in the list, and execute the event with 
earlier occurrence time. We propose the modified version of Direct Method as follows. 
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Algorithm 3. Let {R\,...,Rm} be a set of rewrite rules, X\,...,Xff be numbers of N categories of 
organisms, E be a list of events and maxtime be the time limit for the duration of the simulation. 
Step Initialise simulation time t to 0. Compute propensity a, for every rewrite rule /?,-. 

Step 1 Compute the time increment T. Let (A^e,Ve,?e) be the first event from E with Afe the name of the 
event, Vg the value needed to update the system state and ts the occurrence time of the event. 

Step 2 If t£ < t + T then set t to and then call the event handler algorithm to handle the new event and 
return to Step 1. Otherwise increase simulation time t by time increment z. 

Step 3 If t > maxtime then stop. Otherwise select the next rule index ji and the index 6 of the compart- 
ment in which rule will occur. 

Step 4 Execute rule R^ in the compartment with index 6 and update numbers X\,...,Xn of N cate- 
gories of organisms and propensities a,- for all rewrite rules /?, affected by the application of R^ 
accordingly. Return to Step 1. 

The event handler algorithm is specific to the external events occurring in the system. This algorithm 
updates system state and list of external events and recomputes the propensities that have been affected 
by the change of system state. 

The simulation is affected by the propensity of every rewrite rule. Propensity depends on the number 
of individuals in the population and the rule rate constant. External factors from the environment affect 
propensity values. In general, we cannot associate a rule rate constant with each rewrite rule, because 
the value of the rule rate depends on environmental information, which changes according to external 
events. Since environmental information is incorporated in terms, to model the rule rate we associate 
with that rule a function / ranging over terms. 

Definition 7. Let TV be an infinite set of term variables ranged over by X,Y,Z,. . ., IV be an infinite 
set of information variables ranged over by x,y,z, ■ ■ ■ and NV be an infinite set of natural number vari- 
ables ranged over by q,r,s,... Information Patterns IP, Term Patterns TP and Patterns P are defined as 
follows: 

IP ::= I | I\x CP ::= S \ (T)^ J P TP ::= CP q \ TP\TP P ::= TP \ TP\X 

where X £ TV x £ IV and q £ NV. We denote with £P the set of all patterns. We denote with Var(P) the 
set of variables in P. 

Definition 8. The instantiation is a partial function o : TVUIV UNV -> such that a (TV) C 

g(TV) C ,y and o(NV) C N. We denote with E the set of all possible term instantiations. Given 
P £ &,we denote with Pa the term obtained by replacing all variables X € Var(P) with o(X). 
Definition 9. A rewrite rule is a 4-tuple (f c ,PL,PR,f), usually written as 

[fc]P L ^PR 

where f c : E -> {true, false}, Var(P R ) C Var(P L ), and f :T -> K^°. 

The left pattern matches a portion of the term that models the system by using an instantiation func- 
tion This portion of the system must also satisfy the constraint function f c to enable the rule to be 
applied. A rate function / associated with the rule will be applied to Pr,o. After the rule is applied, P^o 
is substituted by PrO. 

Definition 10. An ecosystem is a triple (T,ffl,E), where T is a term representing the initial state of the 
system, & is a set of rewrite rules that represent the potential internal events which may occur in the 
system, and E is a list of external events. 
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4 Modelling the Population Dynamics of Aedes albopictus 

We use the formalism developed in the previous section to model Aedes albopictus population dynamics. 

4.1 Modelling Information about a Mosquito 

We model each mosquito by using a looping and containment operator with a parallel composition of 
symbols representing information about the mosquito in the content part and a symbol a in the loop part. 
The information in the content part consists of the current development phase of the mosquito and an 
indicator of whether the mosquito has sucked blood or not. In our approach we only model females, 
assuming equal numbers of males and females in the population. In this way we do not need to model 
gender in the information of a mosquito. 

Aedes albopictus goes through 4 development phases in its life cycle: egg, larva, pupa, and adult. 
The larval stage is divided into 4 instars [8]. The adult stage is divided into 8 gonotrophic cycles lfl4"ll . 
A gonotrophic cycle is a cycle in the adult life which consists of three phases called Beklemishev 
phases lfT8l : search for a host and blood-feeding, digestion of the blood and egg maturation, search 
for a suitable oviposition site and oviposition. We use symbols Egg, Larva, Pupa and Adult to denote 
the 4 development phases. Since larva phase is divided into 4 instars, we use symbols 1, 4 to represent 
instars. Analogously, we use symbols 1, 8 to represent gonotrophic cycles. 

An adult mosquito needs blood before ovipositing eggs. We model this phenomenon by adding 
symbol Blood to the content part of the looping and containment operator defining the mosquito to 
represent an adult mosquito that has sucked blood. The number of Blood symbols in the content part 
indicates how many times that mosquito has sucked blood. For instance we represent 3 adult mosquitoes 
at gonotrophic cycle 1 that have sucked blood twice and 5 larvae at instar 1 phase using the following 
term: 

{{a)\ J {Adult | 1| Blood 2 )) 3 | {{a)\ J {Larva | l)) 5 . 

4.2 Modelling Compartments 

In Stochastic CLS compartments are modelled by using looping-containment operators. As we have seen 
in Definition [6] compartments play an important role in our approach, because environmental information 
is attached to them. 

Aedes albopictus, like other species of mosquitoes, spends its immature stages in water. In particular, 
Aedes albopictus prefers to lay eggs outdoors ifTTl . Its natural breeding places are small, restricted, and 
shaded water collections surrounded by vegetation. In urban areas, many man-made containers such as 
tin cans, pots, tires and bottles are usually stored outdoors and collect rainfall water, and thus become 
ideal breeding places lTT5l . Adult Aedes albopictus needs to suck blood before ovipositing. However, 
Aedes albopictus only sucks blood during daytime. Moreover, during immature stages, the duration of 
the stage is affected by temperature while death rate is affected by population density. We can therefore 
define an outermost environmental compartment (we call it environment), with the value of average 
daily temperature and daytime/nighttime as relevant environmental information, inside which there are 
several other compartments where immature mosquitoes live (we call them containers). Population 
density in one container is defined as the number of individuals inside the container divided by the 
water volume in the container. Therefore, relevant environmental information for a container includes 
not only temperature but also water volume and desiccation time. Typical external events are sunrise 
and sunset, which determine switching between daytime and nighttime, temperature changes, which 
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affect desiccation time by reducing the volume of water inside the containers and, as a result, increases 
population density, and rainfalls, which increase the level of water in containers where mosquitoes live, 
so decreasing the population density 

Each kind of compartment has different environmental information. The outermost compartment 
is the environment, to which we need to attach information about current temperature and daylight. 
Therefore, environment is modelled by a term 

(Eft)T em p:V Temp Daylight :V DayUgh , J (T) 

where Vjemp is a real number representing the current temperature, VoayUght is a boolean representing 
whether it is daylight time and T is the term representing the population of Aedes albopictus. 

Immature Aedes albopictus live in small containers, modelled by using looping-containment oper- 
ators with symbol C inside the loop part. For each container we attach the following environmental 
information: 

• an index to identify each container, to be used for container selection by Algorithm |3j 

• the volume of water inside the container, to be used to compute population density; 

• container temperature; 

• three population density thresholds, to be used in the computation of death rates of mosquitoes 
living in the container; 

• container desiccation time. 

If Nc is the number of containers in our model, we use natural numbers in [l,JVc] to identify containers. 
We model the volume of water in an abstract way by classifying containers as full, half —full and 
empty. Population density thresholds, which are used to classify the population density in a container 
and set the death rates accordingly will be defined in Section 14.31 Desiccation, or decrease of water, in 
a container is a process that depends on the characteristic of the container. A desiccation time, which 
measures how many days are needed to reduce the volume of water in a container, is assigned to each 
container. Container desiccation time will be defined in Section |4~4l As an example, term 

Containers ::= (C)f nd . A Temp - W Vohempty fa-.lOO 2 :25O 3 :3OO DTime:2.0 J e I 
(C)fnd:2 TempAO Vol: full 0,:5O <fe:125 3 :15O DTimeA.O J £ 

defines two containers, one identified by number 1, with no water, population density thresholds 100, 
250 and 300, and desiccation time 2 days, and one identified by number 2, full of water, with population 
density thresholds 50, 125 and 150, and desiccation time 1 day. 

A population of immature and adult Aedes albopictus individuals is modelled as a parallel composi- 
tion of looping and containment operators, each with symbol C inside the loop part to model a specific 
container and a parallel composition of looping and containment operators (with symbol a inside the loop 
part) inside the content part to model the immature mosquitoes living inside that container, and looping 
and containment operators with symbol a inside the loop part to model the adult Aedes albopictus in- 
dividuals living in open space. The whole population is then put inside another looping-containment 
operator with symbol En inside the loop part, which models the environment in which the population 
lives. In this way we model the environment in which a population lives as the outermost compartment 
of the Stochastic CLS term that models the biological system of interest. 
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Given the two containers defined above, a daytime environment at a temperature of 10° C with a 
population of 8 adult mosquitoes at the first gonotrophic cycle, 5 of which have sucked blood twice and 
3 of which haven't sucked blood, and 2 empty containers is defined as follows. 

Pop ::= (En)j emp . WDayIjght:tnie J (AdultPop \ Containers) 
AdultPop ::= {{a) L x J (Adult\l\Blood 2 )) 5 \ {a)\ J (Adult\l)f 

We assume that the temperature in all containers is the same as the temperature in the environment. Prop- 
agations of temperature changes in the environment to the containers are handled by the event handler 
algorithm as we will explain in Section l4~4l 

4.3 Modelling Internal Events 

We have seen in Section l4~T1 that the lifecycle of Aedes albopictus consists of the following 14 stages: egg, 
larva (instar 1^1), pupa and adult (8 gonotrophic cycles). Internal events describe transitions between 
some of these stages as well as other events occurring at a specific stage. We identify 29 internal events 
and we model the effect of each of them on the system by a rewrite rule: 

Rule Rl egg hatch 

Rules R2-R4 transitions between instars 

Rule R5 pupation 

Rule R6 adult emergence 

Rule R7 blood sucking 

Rules R8-R15 oviposition at each gonotrophic cycle 

Rules R16-29 death at each stage of the life cycle (14 events) 

Rules R1-R5, which model transitions between immature development stages, rule R6, which mod- 
els the transition from the last immature development stage to the first adult stage, and rules R16-R21, 
which model the death events in such stages, are shown in Figure Q] 

The duration of an immature stage depends on temperature and is measured in degree-days. Degree- 
days for each immature stage is defined as the number of days it takes for an individual in that stage to 
develop at 1°C above the minimum temperature for development (MTD) QJ. Following this definition, 
we define the values of temperature in the environmental information as the difference between the actual 
temperature and MTD. 

If di is the average duration of the i-th development stage, then the rate constant of the rule modelling 
the transition from stage i to the next stage is l/dj. This is true if there are no other events occurring 
during this stage. For every immature development stage we define one rule for the transition to the next 
stage and another one for the death event. Rate functions for transitions in immature stages and death 
events are computed by multiplying 1 / di by survivability rate at i-th stage and by death rate at i-th stage, 
respectively. We assume that the sum of death rate and survivability rate at one development stage is 
equal to one. Since the duration of an immature stage depends on temperature, the rate is then multiplied 
by the difference between the current temperature and MTD. The death rate at an immature stage is 
defined locally for each container and depends on the population density of the container. We classify 
the population density in a container into 4 classes of density: sparse, normal, crowded and overcrowded. 
We define 3 thresholds to be used to classify density: (j>\ , fa and ^3. These three thresholds are part of the 
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Figure 1 : Rewrite rules for the immature stages of Aedes albopictus 



environmental information attached to each container. The rate functions for rules R1-R6 and R16-R21 
are computed as follows: 

y z>«f,- iTnv v v v \ " " (4) 

DD(i— 15) UIfc[10,Zlj 

where 

• j is the index of the rewrite rule, 

•I = ind:k Temp\Vj emp Vol:Vv i 0i fa'-V^ 03 DTime '.Vmime is the environmental 
information attached to the container to which rule /?/ is applied, 

• VVemp is the container temperature, 

• D/?(y,n,Vy o /,V0 1 ,V0 2 ,V^ 3 ) is the death rate function at immature stage j for the container which 
contains n immature mosquitoes, with density thresholds 0i , 02 , 03 and contains a volume Vy / of 
water, 

• DD{j) represents the duration of stage j in degree-days. 

We use 4 classes of population density (sparse, normal, crowded and overcrowded) to define death 
rate in our model. We use the following assumptions for all containers: 

• threshold values used to classify population density in a container are defined for the case in which 
the container is full of water, 

• the baseline death rate of stage i is the death rate of the population in a container whose population 
density is normal, 



28 



Modelling the Dynamics of an Aedes albopictus Population 



• when population in one container is overcrowded or there is no more water in the container only 
death events can occur, so the death rate is set to 1, 

• death rate increases by 20% above the baseline death rate if population density is crowded, 

• death rate decreases by 20% below the baseline death rate if population density is sparse, 

• when a container is only half full, the values of thresholds used to classify the population density 
are divided by 2. 

We define the death rate function D/?:NxNx<fxNxNxN-»Ras follows: 



£>fl(./>,V,01,02,03) = < 



1 
1 

1.2 -BDR(j) 
BDR(j) 
0.&BDR(j) 
, D/?(j',2n,/«//,0i,02,03) 



if V is empty 

if V is full and n > 03 

if V is full and 02 < « < 03 

if V is full and 0i < « < 02 

if V is full and n < 

if V is half —full 



(5) 



where BDR(j) is the baseline death rate for phase j of the life cycle, n is the number of immature 
mosquitoes in the container, 01,02,03 are the container density thresholds and V is the volume of water 
in the container. 

The adult life of an Aedes albopictus is divided into 8 gonotrophic cycles. Every gonotrophic cy- 
cle consists of two internal events: blood sucking and oviposition. The oviposition (the sixth internal 
event) is also a transition from one gonotrophic cycle to the next gonotrophic cycle. Figure [2] shows 
the rewrite rules modelling blood-sucking and oviposition events. Rule Rl models the blood sucking by 
adult mosquitoes. All adult mosquitoes have the same probability of sucking blood. We assume that a 
mosquito always sucks a constant amount of blood. To oviposit, the amount of blood sucked by an adult 
female must be above a threshold (represented by (p in rules R8-R15). 

Rules R8-R15 model the oviposition for the 8 gonotrophic cycles of the mosquito. We assume that 
all adults die after ovipositing at the 8th gonotrophic cycle. The number of eggs any female can produce 
in each gonotrophic cycle is between 45 and 80. This number declines over age. We model this by 
defining function eggs, for each gonotrophic cycle j of the mosquito. 



eggs(j) = < 



' 40 


if 7 


= l 


30 


if j = 5 


37 


if j 


= 2 


27 


if 7 = 6 


35 


if 7 


= 3 


25 


if 7 = 7 


> 32 


if; 


= 4 


22 


if 7 = 8 



(6) 



Although the number of eggs produced by a female mosquito at the j'-th gonotrophic cycle is between 45 
and 80, eggs(j) only returns half of this value to take into account that we only model female individuals. 

Finally figure [3] shows rules R22-R29, which model the death at every adult stage. Rule rates for 
rules R7-R15 and R22-R29 are defined as follows: 



( J_ 



fi 



if i = 7 



I™«2 if ^[8, 15] 



(7) 



where d(i) is the duration of stage i and BDR(i) is the death rate at stage i. 

All rules presented in this section are implemented by using Maude rewrite laws. 
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Figure 2: Rewrite rules for blood-sucking and oviposition events of Aedes albopictus 



4.3.1 Implementation Strategies 

Since we use Maude to implement our model, which mosquito is chosen in the application of rule R7 
depends on the strategy implemented in Maude. To guarantee fairness we implement our own strategy 
in choosing the mosquito with the smallest number of blood sucking times first. 

All adult mosquitoes in a given development stage that have sucked enough blood have the same 
probability of ovipositing. Therefore we consider one rule for each development stage (rules R8-R15). 
We have to deal with the same problem (of choosing the mosquito to oviposit) as in rule R7. To guarantee 
fairness we define a strategy of choosing the mosquito, based on how many times the mosquito has 
sucked blood. As a consequence of the strategy defined for rule R7 the number of times a mosquito 
sucks blood is proportional to the time spent in adult stages. Our strategy will choose the mosquito with 
the biggest number of blood sucking times of ovipositing first. 

We also implement a strategy for choosing the container in which a mosquito oviposits. This strategy 
randomly chooses the container in which the mosquito oviposits. 

The three strategies we have defined in this section have a different purpose from the strategy defined 
by Basuki, Cerone and Milazzo [7], which was used to choose which rewrite rule to apply during a 
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Figure 3: Rewrite rules for death events in adult phases of Aedes albopictus life cycle 

simulation. Instead, the strategies defined in this section are used to choose which portion of the term 
that models the system state matches the lefthand side of a rewrite rule. 

4.4 Modelling External Events 

External events are events that cannot be controlled by the system. These events are usually used to 
model changes in the environment that affect the population. Every event is modelled as a triplet (N,Vt), 
where the event name N is used to distinguish the kind of event, the event value V is used to update the 
environmental information in the system state and the event time t is the time when the event is scheduled 
to occur. Event names and values will be explained in the next paragraphs. Event time t is a non-negative 
real number and measures time in days. The integer part of t represents the day and the fractional part 
represents the time of the day at which an event should occur. For instance t = 1.5 means that the event 
is scheduled to occur on day 1 at 12 pm, and t = 4.125 means that the event is scheduled to occur on day 
4 at 3 am. 

As explained in Section 14.21 for each container there are seven kinds of environmental information 
in our model: container index, container temperature, volume of water in the container, three container 
thresholds for population density and container desiccation time. External events must deal with these 
kinds of environmental information. We define four kinds of event: light change event, change of tem- 
perature, desiccation, and rainfall. Light change events are scheduled twice a day, one at sunrise and 
another at sunset. 

A sunrise event changes the Daylight information associated with the environment from false to 
true. A sunset event changes the Daylight information from true to false. A change of temperature 
event updates temperature in all compartments. A desiccation event updates the volume of water in a 
specific container. A rainfall event updates the volume of water in all containers. Container indices are 
used by the event handler algorithm to handle all events that occur. Population density thresholds are 
used to compute propensity after population density in one container is updated due to the occurrence of 
a desiccation or a rainfall event. Container desiccation time is used to schedule new desiccation events 
due to the occurrence of a desiccation or a rainfall event. 

A light change event is modelled as a triplet {Light, V,t). The time when the sun rises and the time 
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when the sun sets depend on the position of a place on the earth and the time of the year. Value V 
determines whether the event is sunrise (V = sunrise) or sunset (V = sunset) event. For instance in a 
place where in a winter day the sun rises at 8 am and sets at 5 pm, the sunrise event on day 1 is modelled 
as a triplet (Light, sunrise, 1.33) and the sunset event on the same day is modelled as (Light, sunset, 1. .71). 

Temperature affects the duration of immature phases of the mosquito development. We model a tem- 
perature change as a triplet (Temp,Vj emp ,t), which is interpreted as the event of setting the temperature 
to a new value Vjemp starting from time t. We consider only the average daily temperature. We schedule 
one temperature change event every day at midnight. So a triplet (Temp, 10,3.0) means that the average 
temperature on day 3 is 10°C above the MTD of Aedes albopictus. 

The desiccation event is modelled as a triplet (Desic,i,t) which is interpreted as a desiccation in the 
container with index i at time t. We assume that the desiccation time depends on container type and 
measure this time as the number of days needed to reduce the water volume by one level (from full to 
half— full or from half— full to empty). Initially, we introduce one event for each container in list L 
scheduled according to the desiccation time of the container to which it refers. Every time a desiccation 
event occurs and the container is not yet empty, another desiccation event is scheduled to reduce the 
water volume to the next level. For instance, if the system state is represented as: 

(En)f J ((C)f nd . A vohempty DTime:2.0 V J T' \(C)tid:2 Vohfull DTime:\.5 I" J T") 

where / and /" represent part of environmental information which is not relevant for desiccation events, 
T', T" are terms representing population of Aedes albopictus inside container 1 and 2, respectively, and 
the first event in list E is (Desic, 2, 1.0), then at time 1.0 the system state becomes 

(En)f J ((C)f ndA Vohempty DTime:2.Q I' J \ (C)f nd:2 Vohhalf-full DTimeA.S I" J ■ 

The event (Desic, 2, 1.0) is removed from and a new desiccation event (Desic, 2,2.5) is added to list E. 

In our model we only consider containers stored outdoors. In this way, rainfalls are scheduled events 
that increase the water volume level in all containers. Rainfalls are assumed to be prescheduled initially. 
Every time a rainfall event occurs, all desiccation events have to be removed from the list and new 
desiccation events should be added. We classify rainfalls as heavy and light. A heavy rainfall increases 
the water volume level of all containers to full. A light rainfall increases the water volume level of all 
containers from empty to half —full or from half —full to full. The rainfall event is modelled as a 
triplet (Rain, lev,t) which represents a rainfall event with level lev (heavy or light) starting at time t. For 
instance, if the system state is represented as: 

(En)f J ((C)f ndA Vohempty DTime:2.Q I' J ^ \(C)fnd:2 Vohhalf-full DTimeA.5 I" J ^ ) 

and list E contains three events (Rain, light ,1.25), (Desic, 2, 1.5) and (Desic, 1,2.0), then at time 1.25 
the system state becomes 

(En)f J ((C)f ndA Vohhalf-full DTime:2.Q I' J \ (C)f nd:2 Vohfull DTime:\.5 I" J 

The three events are removed from the list and two new desiccation events (Desic, 2,2.1 5) and 
(Desic, 1,3.25) are added to the list. 

The event handler algorithm is very simple. Given a list E of events, Nq containers and a term T that 
represents the system state, the algorithm removes the first event (Ne^e^e) from E and performs the 
different actions decribed above according to the value of Ne £ {Light, Temp, Desic, Rain}. The removal 
of the first event from the list and the subsequent actions are implemented by using Maude rewrite laws. 
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Figure 4: Temperature and Rainfall in Massa Carrara, Italy 



4.5 In silico Experiment and Analysis 

As already mentioned, we have implemented our model in Maude. We have then run a simulation by 
using data collected during May-November 2009 in the province of Massa-Carrara (Tuscany, Italy) in 
1 1 CO2 mosquito traps. The 1 1 traps have captured a total of 3535 Aedes albopictus individuals, and 
have been checked at the following dates: 8 May (4 Aedes a.), 15 May (25 Aedes a.), 19 May (81 Aedes 
a.), 5 June (33 Aedes a.), 18 June (167 Aedes a.), 3 July (360 Aedes a.), 14 July (561 Aedes a.), 29 July 
(381 Aedes a.), 19 August (486 Aedes a.), 3 September (471 Aedes a.), 19 September (276 Aedes a.), 
23 September (292 Aedes a.), 14 October (398 Aedes a.). Note that traps need to be charged with CO2 
in order to work, and that the charge allows the trap to work for one day. Hence, data refer to captures 
of mosquitoes in one day for each considered date. This way of sampling mosquito populations follows 
standard practice. 

Figure 0] shows the climatic data (temperature in °C and rainfalls in mm) during May-November 
2009 in Massa-Carrara province. In our simulation we use 8.8°C as MTD ll24l and 11 containers. Each 
container has carrying capacity of 100-250 organisms and desiccation time between 4.5 and 9.0 days. 

In our simulation we initialise the population with 4 adult mosquitoes (which equals the number of 
adult mosquitoes collected on 8 May 2009) and 10 immature mosquitoes in each of the 11 containers, 
6 eggs, 2 instar- 1 larvae, 1 instar-2 larva and 1 instar-3 larva. The water volume level in each container 
is initially set to half-full. We also set initial desiccation events according the desiccation times of the 
containers. Let to be the time when the simulation starts and D7] be the desiccation time of container 
identified by index i, then we set initial desiccation events at time to + DTj for i = 1 to 1 1 . 

Figure [5] shows the result of our simulation compared with the population sampling produced by us- 
ing the 1 1 traps. We can notice some differences between the simulation results and the field sampling. 
For example, the number of mosquitoes in the sampling decreases between 19 May and 5 June, whereas 
in the simulation such number rapidly increases. This probably happens because of the coarse classifi- 
cation of rainfalls in our model: a very tiny rainfall, with neglectable effect in reality, which occurs just 
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Figure 5: Comparison of in silico simulation (dark line) with data sampled from mosquito traps (light 
line) 



before 19 May, is classified as light rain and, as a result, increases the level of water of the containers in 
the simulation. This may indicate that we need to improve our model by using a finer classification of 
rainfalls. 

The number of mosquitoes captured in traps rapidly increases from 18 June to 14 July, probably due 
to rainfalls. However, no population growth is shown by the simulation during that period. This may be 
due to an overweighed effect that temperature decrease has in our model on immature stage duration and 
death rates. It may also be due to too small values for dessication times used in our model. 

In the simulation the effect of the heavy rainfalls that occur just before 19 September immediately 
causes a population increase on 23 September, but the subsequent decrease in rainfalls and heavy de- 
crease in temperature lead to a population decrease on 14 October. In the field sampling, instead, pop- 
ulation samples keep increasing from 19 September to 14 October. This difference might indicate that 
decrease of rainfalls and temperature take a longer time in reality to affect the population growth than 
in our simulation. This might be again due to too small values for dessication times used in our model. 
Moreover, a decrease in temperature might cause a slower desiccation, a phenomenon that is not consid- 
ered in our model. 



5 Conclusions 

We have presented an extension of the Calculus of Looping Sequences aimed at describing population 
dynamics and ecosystems. The extension consists in allowing a list of external events to be provided 
by the modeller in order to describe environmental events such as changes in the climatic conditions. 
The modeller has also to provide an event handler algorithm that is used by the simulation algorithm 
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associated with the extended formalism. The event handler algorithm is invoked every time an external 
event is planned to occur and it changes the simulation state in accordance with the type of the considered 
event. 

We have used the extended formalism to give a model of a population of Aedes albopictus, an ag- 
gressive mosquito that is well known as an important disease vector. A simulator for this model has been 
developed via translation into the specification language Maude. We have compared results of simula- 
tions of our model with real data obtained from the sampling of mosquitoes during May-November 2009 
in the province of Massa-Carrara (Tuscany, Italy). Since changes in the temperature and rainfalls have a 
significant effect on the mosquito population dynamics, we have exploited data on such environmental 
events (in the same area and the same period of the sampling of mosquitoes) to construct a list of external 
events for the model. 

The results of our simulations show some differences from the real data. However, these differences 
seem to be motivated by some restrictive modelling choices that could be revised in order to construct an 
improved and finer model. Improvements to the model are hence part of our future work, which includes 
also 

• modelling of populations of other disease vector mosquitos such as Aedes aegypti; 

• study of dynamics of populations in other geographic areas; 

• study of different control policies to the mosquito population. 

It would be particularly interesting to study the effects of events such as periodic cleaning of containers 
and use of pesticides on the mosquito population to choose the most promising control policy. Such a 
policy could then be experimented in the field, and the results obtained could be used to further validate 
the model. A method to choose the best mosquito population control policy would be of interest in 
particular in those areas in which such mosquitoes act as vectors of diseases. 
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